Mechanistic modeling of climate effects on redistribution and population growth in a community of fish species.

Understanding community responses to climate is critical for anticipating the future impacts of global change. However, despite increased research efforts in this field, models that explicitly include important biological mechanisms are lacking. Quantifying the potential impacts of climate change on species is complicated by the fact that the effects of climate variation may manifest at several points in the biological process. To this end, we extend a dynamic mechanistic model that combines population dynamics, such as species interactions, with species redistribution by allowing climate to affect both processes. We examine their relative contributions in an application to the changing biomass of a community of eight species in the Gulf of Maine using over 30 years of fisheries data from the Northeast Fishery Science Center. Our model suggests that the mechanisms driving biomass trends vary across space, time, and species. Phase space plots demonstrate that failing to account for the dynamic nature of the environmental and biologic system can yield theoretical estimates of population abundances that are not observed in empirical data. The stock assessments used by fisheries managers to set fishing targets and allocate quotas often ignore environmental effects. At the same time, research examining the effects of climate change on fish has largely focused on redistribution. Frameworks that combine multiple biological reactions to climate change are particularly necessary for marine researchers. This work is just one approach to modeling the complexity of natural systems and highlights the need to incorporate multiple and possibly interacting biological processes in future models.

Assessing Disparity Using Measures of Racial and Educational Isolation.

We develop a local, spatial measure of educational isolation (EI) and characterize the relationship between EI and our previously developed measure of racial isolation (RI). EI measures the extent to which non-college educated individuals are exposed primarily to other non-college educated individuals. To characterize how the RI-EI relationship varies across space, we propose a novel measure of local correlation. Using birth records from the State of Michigan (2005-2012), we estimate associations between RI, EI, and birth outcomes. EI was lower in urban communities and higher in rural communities, while RI was highest in urban areas and parts of the southeastern United States (US). We observed greater heterogeneity in EI in low RI tracts, especially in non-urban tracts; residents of high RI tracts are likely to be both educationally and racially isolated. Associations were also observed between RI, EI, and gestational length (weeks) and preterm birth (PTB). For example, moving from the lowest to the highest quintile of RI was associated with a 1.11 (1.07, 1.15) and 1.16 (1.10, 1.22) increase in odds of PTB among NHB and NHW women, respectively. Moving from the lowest to the highest quintile of EI was associated with a 1.07 (1.02, 1.12) and 1.03 (1.00, 1.05) increase in odds of PTB among NHB and NHW women, respectively. This work provides three tools (RI, EI, and the local correlation measure) to researchers and policymakers interested in how residential isolation shapes disparate outcomes.

Spatial Joint Species Distribution Modeling using Dirichlet Processes.

Species distribution models usually attempt to explain presence-absence or abundance of a species at a site in terms of the environmental features (so-called abiotic features) present at the site. Historically, such models have considered species individually. However, it is well-established that species interact to influence presence-absence and abundance (envisioned as biotic factors). As a result, there has been substantial recent interest in joint species distribution models with various types of response, e.g., presence-absence, continuous and ordinal data. Such models incorporate dependence between species response as a surrogate for interaction. The challenge we address here is how to accommodate such modeling in the context of a large number of species (e.g., order 102) across sites numbering on the order of 102 or 103 when, in practice, only a few species are found at any observed site. Again, there is some recent literature to address this; we adopt a dimension reduction approach. The novel wrinkle we add here is spatial dependence. That is, we have a collection of sites over a relatively small spatial region so it is anticipated that species distribution at a given site would be similar to that at a nearby site. Specifically, we handle dimension reduction through Dirichlet processes, enabling clustering of species, joined with spatial dependence across sites through Gaussian processes. We use both simulated data and a plant communities dataset for the Cape Floristic Region (CFR) of South Africa to demonstrate our approach. The latter consists of presence-absence measurements for 639 tree species at 662 locations. Through both data examples we are able to demonstrate improved predictive performance using the foregoing specification.

Alternating Gaussian Process Modulated Renewal Processes for Modeling Threshold Exceedances and Durations.

It is often of interest to model the incidence and duration of threshold exceedance events for an environmental variable over a set of monitoring locations. Such data arrive over continuous time and can be considered as observations of a two-state process yielding, sequentially, a length of time in the below threshold state followed by a length of time in the above threshold state, then returning to the below threshold state, etc. We have a two-state continuous time Markov process, often referred to as an alternating renewal process. The process is observed over a truncated time window and, within this window, time in each state is modeled using a distinct cumulative intensity specification. Initially, we model each intensity over the window using a parametric regression specification. We extend the regression specification adding temporal random effects to enrich the model, using a realization of a log Gaussian process over time. With only one type of renewal, this specification is referred to as a Gaussian process modulated renewal process. Here, we introduce Gaussian process modulation to the intensity for each state. Model fitting is done within a Bayesian framework. We clarify that fitting with a customary log Gaussian process specification over a lengthy time window is computationally infeasible. The nearest neighbor Gaussian process (NNGP), which supplies sparse covariance structure, is adopted to enable tractable computation. We also propose methods for both generating data under our models and for conducting model comparison. The model is applied to hourly ozone data for four monitoring sites in different locations across the United States for the ozone season of 2014. For each site, we obtain estimated profiles of up-crossing and down-crossing intensity functions through time. In addition, we obtain inference regarding the number of exceedances, the distribution of the duration of exceedance events, and the proportion of time in the above and below threshold state for any time interval.

Accommodating the ecological fallacy in disease mapping in the absence of individual exposures.

In health exposure modeling, in particular, disease mapping, the ecological fallacy arises because the relationship between aggregated disease incidence on areal units and average exposure on those units differs from the relationship between the event of individual incidence and the associated individual exposure. This article presents a novel modeling approach to address the ecological fallacy in the least informative data setting. We assume the known population at risk with an observed incidence for a collection of areal units and, separately, environmental exposure recorded during the period of incidence at a collection of monitoring stations. We do not assume any partial individual level information or random allocation of individuals to observed exposures. We specify a conceptual incidence surface over the study region as a function of an exposure surface resulting in a stochastic integral of the block average disease incidence. The true block level incidence is an unavailable Monte Carlo integration for this stochastic integral. We propose an alternative manageable Monte Carlo integration for the integral. Modeling in this setting is immediately hierarchical, and we fit our model within a Bayesian framework. To alleviate the resulting computational burden, we offer 2 strategies for efficient model fitting: one is through modularization, the other is through sparse or dimension-reduced Gaussian processes. We illustrate the performance of our model with simulations based on a heat-related mortality dataset in Ohio and then analyze associated real data.

Bayesian Modeling and Analysis of Geostatistical Data.

The most prevalent spatial data setting is, arguably, that of so-called geostatistical data, data that arise as random variables observed at fixed spatial locations. Collection of such data in space and in time has grown enormously in the past two decades. With it has grown a substantial array of methods to analyze such data. Here, we attempt a review of a fully model-based perspective for such data analysis, the approach of hierarchical modeling fitted within a Bayesian framework. The benefit, as with hierarchical Bayesian modeling in general, is full and exact inference, with proper assessment of uncertainty. Geostatistical modeling includes univariate and multivariate data collection at sites, continuous and categorical data at sites, static and dynamic data at sites, and datasets over very large numbers of sites and long periods of time. Within the hierarchical modeling framework, we offer a review of the current state of the art in these settings.

On nearest-neighbor Gaussian process models for massive spatial data.

Gaussian Process (GP) models provide a very flexible nonparametric approach to modeling location-and-time indexed datasets. However, the storage and computational requirements for GP models are infeasible for large spatial datasets. Nearest Neighbor Gaussian Processes (Datta A, Banerjee S, Finley AO, Gelfand AE. Hierarchical nearest-neighbor gaussian process models for large geostatistical datasets. J Am Stat Assoc 2016., JASA) provide a scalable alternative by using local information from few nearest neighbors. Scalability is achieved by using the neighbor sets in a conditional specification of the model. We show how this is equivalent to sparse modeling of Cholesky factors of large covariance matrices. We also discuss a general approach to construct scalable Gaussian Processes using sparse local kriging. We present a multivariate data analysis which demonstrates how the nearest neighbor approach yields inference indistinguishable from the full rank GP despite being several times faster. Finally, we also propose a variant of the NNGP model for automating the selection of the neighbor set size.

Hierarchical Nearest-Neighbor Gaussian Process Models for Large Geostatistical Datasets.

Spatial process models for analyzing geostatistical data entail computations that become prohibitive as the number of spatial locations become large. This article develops a class of highly scalable nearest-neighbor Gaussian process (NNGP) models to provide fully model-based inference for large geostatistical datasets. We establish that the NNGP is a well-defined spatial process providing legitimate finite-dimensional Gaussian densities with sparse precision matrices. We embed the NNGP as a sparsity-inducing prior within a rich hierarchical modeling framework and outline how computationally efficient Markov chain Monte Carlo (MCMC) algorithms can be executed without storing or decomposing large matrices. The floating point operations (flops) per iteration of this algorithm is linear in the number of spatial locations, thereby rendering substantial scalability. We illustrate the computational and inferential benefits of the NNGP over competing methods using simulation studies and also analyze forest biomass from a massive U.S. Forest Inventory dataset at a scale that precludes alternative dimension-reducing methods. Supplementary materials for this article are available online.

More than the sum of the parts: forest climate response from joint species distribution models.

The perceived threat of climate change is often evaluated from species distribution models that are fitted to many species independently and then added together. This approach ignores the fact that species are jointly distributed and limit one another. Species respond to the same underlying climatic variables, and the abundance of any one species can be constrained by competition; a large increase in one is inevitably linked to declines of others. Omitting this basic relationship explains why responses modeled independently do not agree with the species richness or basal areas of actual forests. We introduce a joint species distribution modeling approach (JSDM), which is unique in three ways, and apply it to forests of eastern North America. First, it accommodates the joint distribution of species. Second, this joint distribution includes both abundance and presence-absence data. We solve the common issue of large numbers of zeros in abundance data by accommodating zeros in both stem counts and basal area data, i.e., a new approach to zero inflation. Finally, inverse prediction can be applied to the joint distribution of predictions to integrate the role of climate risks across all species and identify geographic areas where communities will change most (in terms of changes in abundance) with climate change. Application to forests in the eastern United States shows that climate can have greatest impact in the Northeast, due to temperature, and in the Upper Midwest, due to temperature and precipitation. Thus, these are the regions experiencing the fastest warming and are also identified as most responsive at this scale.

Imputation of confidential data sets with spatial locations using disease mapping models.

Data that include fine geographic information, such as census tract or street block identifiers, can be difficult to release as public use files. Fine geography provides information that ill-intentioned data users can use to identify individuals. We propose to release data with simulated geographies, so as to enable spatial analyses while reducing disclosure risks. We fit disease mapping models that predict areal-level counts from attributes in the file and sample new locations based on the estimated models. We illustrate this approach using data on causes of death in North Carolina, including evaluations of the disclosure risks and analytic validity that can result from releasing synthetic geographies.

A multivariate spatial mixture model for areal data: examining regional differences in standardized test scores.

Researchers in the health and social sciences often wish to examine joint spatial patterns for two or more related outcomes. Examples include infant birth weight and gestational length, psychosocial and behavioral indices, and educational test scores from different cognitive domains. We propose a multivariate spatial mixture model for the joint analysis of continuous individual-level outcomes that are referenced to areal units. The responses are modeled as a finite mixture of multivariate normals, which accommodates a wide range of marginal response distributions and allows investigators to examine covariate effects within subpopulations of interest. The model has a hierarchical structure built at the individual level (i.e., individuals are nested within areal units), and thus incorporates both individual- and areal-level predictors as well as spatial random effects for each mixture component. Conditional autoregressive (CAR) priors on the random effects provide spatial smoothing and allow the shape of the multivariate distribution to vary flexibly across geographic regions. We adopt a Bayesian modeling approach and develop an efficient Markov chain Monte Carlo model fitting algorithm that relies primarily on closed-form full conditionals. We use the model to explore geographic patterns in end-of-grade math and reading test scores among school-age children in North Carolina.

Dual impacts of climate change: forest migration and turnover through life history.

Tree species are predicted to track future climate by shifting their geographic distributions, but climate-mediated migrations are not apparent in a recent continental-scale analysis. To better understand the mechanisms of a possible migration lag, we analyzed relative recruitment patterns by comparing juvenile and adult tree abundances in climate space. One would expect relative recruitment to be higher in cold and dry climates as a result of tree migration with juveniles located further poleward than adults. Alternatively, relative recruitment could be higher in warm and wet climates as a result of higher tree population turnover with increased temperature and precipitation. Using the USDA Forest Service's Forest Inventory and Analysis data at regional scales, we jointly modeled juvenile and adult abundance distributions for 65 tree species in climate space of the eastern United States. We directly compared the optimal climate conditions for juveniles and adults, identified the climates where each species has high relative recruitment, and synthesized relative recruitment patterns across species. Results suggest that for 77% and 83% of the tree species, juveniles have higher optimal temperature and optimal precipitation, respectively, than adults. Across species, the relative recruitment pattern is dominated by relatively more abundant juveniles than adults in warm and wet climates. These different abundance-climate responses through life history are consistent with faster population turnover and inconsistent with the geographic trend of large-scale tree migration. Taken together, this juvenile-adult analysis suggests that tree species might respond to climate change by having faster turnover as dynamics accelerate with longer growing seasons and higher temperatures, before there is evidence of poleward migration at biogeographic scales.

Modeling daily flowering probabilities: expected impact of climate change on Japanese cherry phenology.

Understanding the drivers of phenological events is vital for forecasting species' responses to climate change. We developed flexible Bayesian survival regression models to assess a 29-year, individual-level time series of flowering phenology from four taxa of Japanese cherry trees (Prunus spachiana, Prunus × yedoensis, Prunus jamasakura, and Prunus lannesiana), from the Tama Forest Cherry Preservation Garden in Hachioji, Japan. Our modeling framework used time-varying (chill and heat units) and time-invariant (slope, aspect, and elevation) factors. We found limited differences among taxa in sensitivity to chill, but earlier flowering taxa, such as P. spachiana, were more sensitive to heat than later flowering taxa, such as P. lannesiana. Using an ensemble of three downscaled regional climate models under the A1B emissions scenario, we projected shifts in flowering timing by 2100. Projections suggest that each taxa will flower about 30 days earlier on average by 2100 with 2-6 days greater uncertainty around the species mean flowering date. Dramatic shifts in the flowering times of cherry trees may have implications for economically important cultural festivals in Japan and East Asia. The survival models used here provide a mechanistic modeling approach and are broadly applicable to any time-to-event phenological data, such as plant leafing, bird arrival time, and insect emergence. The ability to explicitly quantify uncertainty, examine phenological responses on a fine time scale, and incorporate conditions leading up to an event may provide future insight into phenologically driven changes in carbon balance and ecological mismatches of plants and pollinators in natural populations and horticultural crops.

Directional data analysis under the general projected normal distribution.

The projected normal distribution is an under-utilized model for explaining directional data. In particular, the general version provides flexibility, e.g., asymmetry and possible bimodality along with convenient regression specification. Here, we clarify the properties of this general class. We also develop fully Bayesian hierarchical models for analyzing circular data using this class. We show how they can be fit using MCMC methods with suitable latent variables. We show how posterior inference for distributional features such as the angular mean direction and concentration can be implemented as well as how prediction within the regression setting can be handled. With regard to model comparison, we argue for an out-of-sample approach using both a predictive likelihood scoring loss criterion and a cumulative rank probability score criterion.

Spatio-temporal modeling for real-time ozone forecasting.

The accurate assessment of exposure to ambient ozone concentrations is important for informing the public and pollution monitoring agencies about ozone levels that may lead to adverse health effects. High-resolution air quality information can offer significant health benefits by leading to improved environmental decisions. A practical challenge facing the U.S. Environmental Protection Agency (USEPA) is to provide real-time forecasting of current 8-hour average ozone exposure over the entire conterminous United States. Such real-time forecasting is now provided as spatial forecast maps of current 8-hour average ozone defined as the average of the previous four hours, current hour, and predictions for the next three hours. Current 8-hour average patterns are updated hourly throughout the day on the EPA-AIRNow web site. The contribution here is to show how we can substantially improve upon current real-time forecasting systems. To enable such forecasting, we introduce a downscaler fusion model based on first differences of real-time monitoring data and numerical model output. The model has a flexible coefficient structure and uses an efficient computational strategy to fit model parameters. Our hybrid computational strategy blends continuous background updated model fitting with real-time predictions. Model validation analyses show that we are achieving very accurate and precise ozone forecasts.

The k-ZIG: flexible modeling for zero-inflated counts.

Many applications involve count data from a process that yields an excess number of zeros. Zero-inflated count models, in particular, zero-inflated Poisson (ZIP) and zero-inflated negative binomial (ZINB) models, along with Poisson hurdle models, are commonly used to address this problem. However, these models struggle to explain extreme incidence of zeros (say more than 80%), especially to find important covariates. In fact, the ZIP may struggle even when the proportion is not extreme. To redress this problem we propose the class of k-ZIG models. These models allow more flexible modeling of both the zero-inflation and the nonzero counts, allowing interplay between these two components. We develop the properties of this new class of models, including reparameterization to a natural link function. The models are straightforwardly fitted within a Bayesian framework. The methodology is illustrated with simulated data examples as well as a forest seedling dataset obtained from the USDA Forest Service's Forest Inventory and Analysis program.

Space-time data fusion under error in computer model output: an application to modeling air quality.

We provide methods that can be used to obtain more accurate environmental exposure assessment. In particular, we propose two modeling approaches to combine monitoring data at point level with numerical model output at grid cell level, yielding improved prediction of ambient exposure at point level. Extending our earlier downscaler model (Berrocal, V. J., Gelfand, A. E., and Holland, D. M. (2010b). A spatio-temporal downscaler for outputs from numerical models. Journal of Agricultural, Biological and Environmental Statistics 15, 176-197), these new models are intended to address two potential concerns with the model output. One recognizes that there may be useful information in the outputs for grid cells that are neighbors of the one in which the location lies. The second acknowledges potential spatial misalignment between a station and its putatively associated grid cell. The first model is a Gaussian Markov random field smoothed downscaler that relates monitoring station data and computer model output via the introduction of a latent Gaussian Markov random field linked to both sources of data. The second model is a smoothed downscaler with spatially varying random weights defined through a latent Gaussian process and an exponential kernel function, that yields, at each site, a new variable on which the monitoring station data is regressed with a spatial linear model. We applied both methods to daily ozone concentration data for the Eastern US during the summer months of June, July and August 2001, obtaining, respectively, a 5% and a 15% predictive gain in overall predictive mean square error over our earlier downscaler model (Berrocal et al., 2010b). Perhaps more importantly, the predictive gain is greater at hold-out sites that are far from monitoring sites.

Hierarchical Modeling for Spatial Data Problems.

This short paper is centered on hierarchical modeling for problems in spatial and spatio-temporal statistics. It draws its motivation from the interdisciplinary research work of the author in terms of applications in the environmental sciences - ecological processes, environmental exposure, and weather modeling. The paper briefly reviews hierarchical modeling specification, adopting a Bayesian perspective with full inference and associated uncertainty within the specification, while achieving exact inference to avoid what may be uncomfortable asymptotics. It focuses on point-referenced (geo-statistical) and point pattern spatial settings. It looks in some detail at problems involving data fusion, species distributions, and large spatial datasets. It also briefly describes four further examples arising from the author's recent research projects.

Kernel Averaged Predictors for Spatio-Temporal Regression Models.

In applications where covariates and responses are observed across space and time, a common goal is to quantify the effect of a change in the covariates on the response while adequately accounting for the spatio-temporal structure of the observations. The most common approach for building such a model is to confine the relationship between a covariate and response variable to a single spatio-temporal location. However, oftentimes the relationship between the response and predictors may extend across space and time. In other words, the response may be affected by levels of predictors in spatio-temporal proximity to the response location. Here, a flexible modeling framework is proposed to capture such spatial and temporal lagged effects between a predictor and a response. Specifically, kernel functions are used to weight a spatio-temporal covariate surface in a regression model for the response. The kernels are assumed to be parametric and non-stationary with the data informing the parameter values of the kernel. The methodology is illustrated on simulated data as well as a physical data set of ozone concentrations to be explained by temperature.